function IR = LPM(I)
 [F C] = size(I);
 
 Om = (F+1)/2;                  % coordenadas del centro de la imagen
 On = (C+1)/2;
 sx = (F-1)/2;                  % factores de escala
 sy = (C-1)/2;
 
 M=64;
 N=360;
  
 %img=zeros(M*2,N);
 
 R=0.5;
 a=2^(1/M);
 rmin=(a^(-M/2))*R;
 rmax=(a^(M/2))*R;

 cont=zeros(M,N);
 suma=zeros(M,N);
 img=zeros(M,N);
 
for xi = 1:F
  for yi = 1:C
     %u = (xi - Om)/sx;
     %v = (((-1)*yi + On)/sy);
     
     u = ((xi - Om)/sx);
     v = ((yi - On)/sy);
     
     r= sqrt(u^2 + v^2);

     if (r >=rmin) && (r <rmax)
        
         theta = atan2(v,u);
       
        
         l1 =abs(floor((log2(r)-log2(R))/log2(a))+M/2);
         % l1 =abs(floor((log(r/R)/(log(a)))+M/2)); 
         l2 =(floor((N*theta)/pi));
          
        if(l1>=0 && l1<M) && (l2>=0 && l2<N)
          img(l1+1,l2+1) = I(xi+1,yi+1);
           
          cont(l1+1,l2+1)=cont(l1+1,l2+1)+1;
          suma(l1+1,l2+1)=suma(l1+1,l2+1)+I(xi+1,yi+1);

        end
      end 
  end
end
cont(cont(:)==0)=1;
img(:)=suma(:)./cont(:);
IR=img;
end